package flare.physics
{

/**
 * Force simulating an N-Body force of charged particles with pairwise
 * interaction, such as gravity or electrical charge. This class uses a
 * quad-tree structure to aggregate charge values and optimize computation.
 * The force function is a standard inverse-square law (though in this case
 * approximated due to optimization): <code>F = G * m1 * m2 / d^2</code>,
 * where G is a constant (e.g., gravitational constant), m1 and m2 are the
 * masses (charge) of the particles, and d is the distance between them.
 * 
 * <p>The algorithm used is that of J. Barnes and P. Hut, in their research
 * paper <i>A Hierarchical  O(n log n) force calculation algorithm</i>, Nature, 
 * v.324, December 1986. For more details on the algorithm, see one of
 * the following links:
 * <ul>
 *   <li><a href="http://www.cs.berkeley.edu/~demmel/cs267/lecture26/lecture26.html">James Demmel's UC Berkeley lecture notes</a>
 *   <li><a href="http://www.physics.gmu.edu/~large/lr_forces/desc/bh/bhdesc.html">Description of the Barnes-Hut algorithm</a>
 *   <li><a href="http://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html">Joshua Barnes' implementation</a>
 * </ul></p>
 */
  public class NBodyForce implements IForce
  {
    private var _g:Number;   // gravitational constant
    private var _t:Number;   // barnes-hut theta
    private var _max:Number; // max effective distance
    private var _min:Number; // min effective distance
    private var _eps:Number; // epsilon for determining 'same' location
    private var _x1:Number, _y1:Number, _x2:Number, _y2:Number;
    private var _root:QuadTreeNode;

    /** The gravitational constant to use. 
     *  Negative values produce a repulsive force. */
    public function get gravitation():Number
    {
      return _g;
    }

    public function set gravitation(g:Number):void
    {
      _g = g;
    }

    /** The maximum distance over which forces are exerted. 
     *  Any greater distances will be ignored. */
    public function get maxDistance():Number
    {
      return _max;
    }

    public function set maxDistance(d:Number):void
    {
      _max = d;
    }

    /** The minumum effective distance over which forces are exerted.
     * 	Any lesser distances will be treated as the minimum. */
    public function get minDistance():Number
    {
      return _min;
    }

    public function set minDistance(d:Number):void
    {
      _min = d;
    }

    // --------------------------------------------------------------------

    /**
     * Creates a new NBodyForce with given parameters.
     * @param g the gravitational constant to use.
     *  Negative values produce a repulsive force.
     * @param maxd a maximum distance over which the force should operate.
     *  Particles separated by more than this distance will not interact.
     * @param mind the minimum distance over which the force should operate.
     *  Particles closer than this distance will interact as if they were
     *  the minimum distance apart. This helps avoid extreme forces.
     *  Helpful when particles are very close together.
     * @param eps an epsilon values for determining a minimum distance
     *  between particles
     * @param t the theta parameter for the Barnes-Hut approximation.
     *  Determines the level of approximation (default value if 0.9).
     */
    public function NBodyForce(g:Number = - 1, max:Number = 200, min:Number = 2, eps:Number = 0.01, t:Number = 0.9)
    {
      _g = g;
      _max = max;
      _min = min;
      _eps = eps;
      _t = t;
      _root = QuadTreeNode.node();
    }

    /**
     * Applies this force to a simulation.
     * @param sim the Simulation to apply the force to
     */
    public function apply(sim:Simulation):void
    {
      if(_g == 0)return;

      // clear the quadtree
      clear(_root);
      _root = QuadTreeNode.node();

      // get the tree bounds
      bounds(sim);

      // populate the tree
      for(var i:uint = 0; i < sim.particles.length; ++i)
      {
        insert(sim.particles[i], _root, _x1, _y1, _x2, _y2);
      }

      // traverse tree to compute mass
      accumulate(_root);

      // calculate forces on each particle
      for(i = 0; i < sim.particles.length; ++i)
      {
        forces(sim.particles[i], _root, _x1, _y1, _x2, _y2);
      }
    }

    private function accumulate(n:QuadTreeNode):void
    {
      var xc:Number = 0, yc:Number = 0;
      n.mass = 0;

      // accumulate childrens' mass
      var recurse:Function = function(c:QuadTreeNode):void
      {
        if(c == null)return;

        accumulate(c);
        n.mass += c.mass;
        xc += c.mass * c.cx;
        yc += c.mass * c.cy;
      }

      if(n.hasChildren)
      {
        recurse(n.c1);
        recurse(n.c2);
        recurse(n.c3);
        recurse(n.c4);
      }

      // accumulate own mass
      if(n.p != null)
      {
        n.mass += n.p.mass;
        xc += n.p.mass * n.p.x;
        yc += n.p.mass * n.p.y;
      }
      n.cx = xc / n.mass;
      n.cy = yc / n.mass;
    }

    private function forces(p:Particle, n:QuadTreeNode, x1:Number, y1:Number, x2:Number, y2:Number):void
    {
      var f:Number = 0;
      var dx:Number = n.cx - p.x;
      var dy:Number = n.cy - p.y;
      var dd:Number = Math.sqrt(dx * dx + dy * dy);
      var max:Boolean = _max > 0 && dd > _max;

      if(dd == 0)
      {
        // add direction when needed
        dx = _eps * (0.5 - Math.random());
        dy = _eps * (0.5 - Math.random());
      }

      // the Barnes-Hut approximation criteria is if the ratio of the
      // size of the quadtree box to the distance between the point and
      // the box's center of mass is beneath some threshold theta.
      if((!n.hasChildren && n.p != p) || ((x2 - x1) / dd < _t))
      {
        if(max)return;
        // either only 1 particle or we meet criteria
        // for Barnes-Hut approximation, so calc force
        dd = dd < _min ? _min : dd;
        f = _g * p.mass * n.mass / (dd * dd * dd)
        p.fx += f * dx;
        p.fy += f * dy;
      }
      else if(n.hasChildren)
      {
        // recurse for more accurate calculation
        var sx:Number = (x1 + x2) / 2
        var sy:Number = (y1 + y2) / 2;

        if(n.c1)forces(p, n.c1, x1, y1, sx, sy);

        if(n.c2)forces(p, n.c2, sx, y1, x2, sy);

        if(n.c3)forces(p, n.c3, x1, sy, sx, y2);

        if(n.c4)forces(p, n.c4, sx, sy, x2, y2);

        if(max)return;

        if(n.p != null && n.p != p)
        {
          dd = dd < _min ? _min : dd;
          f = _g * p.mass * n.p.mass / (dd * dd * dd);
          p.fx += f * dx;
          p.fy += f * dy;
        }
      }
    }

    // -- Helpers ---------------------------------------------------------
    private function insert(p:Particle, n:QuadTreeNode, x1:Number, y1:Number, x2:Number, y2:Number):void
    {
      // ignore particles with NaN coordinates
      if(isNaN(p.x) || isNaN(p.y))return;

      // try to insert particle p at node n in the quadtree
      // by construction, each leaf will contain either 1 or 0 particles
      if(n.hasChildren)
      {
        // n contains more than 1 particle
        insertHelper(p, n, x1, y1, x2, y2);
      }
      else if(n.p != null)
      {
        // n contains 1 particle
        if(isSameLocation(n.p, p))
        {
          // recurse
          insertHelper(p, n, x1, y1, x2, y2);
        }

        else
        {
          // divide
          var v:Particle = n.p;
          n.p = null;
          insertHelper(v, n, x1, y1, x2, y2);
          insertHelper(p, n, x1, y1, x2, y2);
        }
      }

      else
      {
        // n is empty, add p as leaf
        n.p = p;
      }
    }

    private function insertHelper(p:Particle, n:QuadTreeNode, x1:Number, y1:Number, x2:Number, y2:Number):void
    {
      // determine split
      var sx:Number = (x1 + x2) / 2;
      var sy:Number = (y1 + y2) / 2;
      var c:uint = (p.x >= sx ? 1 : 0) + (p.y >= sy ? 2 : 0);

      // update bounds
      if(c == 1 || c == 3)x1 = sx;

      else x2 = sx;

      if(c > 1)y1 = sy;

      else y2 = sy;

      // update children
      var cn:QuadTreeNode;

      if(c == 0)
      {
        if(n.c1 == null)n.c1 = QuadTreeNode.node();
        cn = n.c1;
      }
      else if(c == 1)
      {
        if(n.c2 == null)n.c2 = QuadTreeNode.node();
        cn = n.c2;
      }
      else if(c == 2)
      {
        if(n.c3 == null)n.c3 = QuadTreeNode.node();
        cn = n.c3;
      }

      else
      {
        if(n.c4 == null)n.c4 = QuadTreeNode.node();
        cn = n.c4;
      }
      n.hasChildren = true;
      insert(p, cn, x1, y1, x2, y2);
    }

    private function clear(n:QuadTreeNode):void
    {
      if(n.c1 != null)clear(n.c1);

      if(n.c2 != null)clear(n.c2);

      if(n.c3 != null)clear(n.c3);

      if(n.c4 != null)clear(n.c4);
      QuadTreeNode.reclaim(n);
    }

    private function bounds(sim:Simulation):void
    {
      var p:Particle, dx:Number, dy:Number;
      _x1 = _y1 = Number.MAX_VALUE;
      _x2 = _y2 = Number.MIN_VALUE;

      // get bounding box
      for(var i:uint = 0; i < sim.particles.length; ++i)
      {
        p = sim.particles[i] as Particle;

        if(p.x < _x1)_x1 = p.x;

        if(p.y < _y1)_y1 = p.y;

        if(p.x > _x2)_x2 = p.x;

        if(p.y > _y2)_y2 = p.y;
      }

      // square the box
      dx = _x2 - _x1;
      dy = _y2 - _y1;

      if(dx > dy)
      {
        _y2 = _y1 + dx;
      }

      else
      {
        _x2 = _x1 + dy;
      }
    }

    private function isSameLocation(p1:Particle, p2:Particle):Boolean
    {
      return (Math.abs(p1.x - p2.x) < _eps && Math.abs(p1.y - p2.y) < _eps);
    }
  } // end of class NBodyForce
}

import flare.physics.Particle;

// -- Helper QuadTreeNode class -----------------------------------------------

class QuadTreeNode
{

  public var mass:Number = 0;
  public var cx:Number = 0;
  public var cy:Number = 0;
  public var p:Particle = null;
  public var c1:QuadTreeNode = null;
  public var c2:QuadTreeNode = null;
  public var c3:QuadTreeNode = null;
  public var c4:QuadTreeNode = null;
  public var hasChildren:Boolean = false;
  // -- Factory ---------------------------------------------------------
  private static var _nodes:Vector.<Object> = new Vector.<Object>();

  public static function node():QuadTreeNode
  {
    var n:QuadTreeNode;

    if(_nodes.length > 0)
    {
      n = QuadTreeNode(_nodes.pop());
    }

    else
    {
      n = new QuadTreeNode();
    }
    return n;
  }

  public static function reclaim(n:QuadTreeNode):void
  {
    n.mass = n.cx = n.cy = 0;
    n.p = null;
    n.hasChildren = false;
    n.c1 = n.c2 = n.c3 = n.c4 = null;
    _nodes.push(n);
  }
}